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PREFACE 


Part  of  the  RAND  research  program  consists  of  basic 
supporting  studies  in  mathematics.  This  Memorandum  is 
the  first  in  a  series  dealing  with  a  number  of  rigorous 
aspects  of  the  highly  useful  mathematical  theory  known 
as  invariant  imbedding.  In  this  theory  invariance  prin¬ 
ciples  are  applied  to  handle  a  variety  of  conceptual  and 
computational  aspects  of  mathematical  physics. 

The  research  presented  here  was  sponsored  by  the 
Advanced  Research  Projects  Agency. 
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SUMMARY 

In  a  series  of  papers  of  which  this  is  the  first, 
we  wish  to  study  some  of  the  rigorous  aspects  of  invariant 
imbedding:  existence  and  uniqueness  of  solution,  asympto¬ 
tic  behavior  over  space  and  time,  stability,  computational 
stability,  applications  to  classical  boundary— value  theory, 
and  so  on. 

The  first  paper  will  be  devoted  to  the  use  of  an 
important  conservation  property,  obvious  on  physical 
grounds,  in  establishing  the  existence  of  the  solution  of 
a  matrix  Riccati  equation  without  recourse  to  the  asso¬ 
ciated  linear  differential  equation,  and  thus  without 
any  appeal  to  spectral  theory. 
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EXISTENCE  AND  UNIQUENESS  THEOREMS  IN  INVARIANT  IMBEDDING  - 
Is  CONSERVATION  PRINCIPLES 


1 .  INTRODUCTION 

Invariant  imbedding  is  a  mathematical  theory  designed 
to  handle  a  variety  of  conceptual,  analytic,  and  computa¬ 
tional  aspects  of  mathematical  physics  in  a  unified  fashion 
without  the  intervention  of  boundary— value  problems.  By 
means  of  appropriate  choices  of  space  and  time  variables, 
all  problems  are  of  initial— value  type.  An  expository 
account  of  the  application  of  this  theory  to  neutron  trans¬ 
port,  radiative  transfer,  diffusion,  and  scattering  will 
be  found  in  [i];  application  to  wave  propagation  will  be 
found  in  [2],  [3]. 

Invariant  Imbedding  is  a  systematic  application  and 
extension  of  the  "invariance  principles"  introduced  into 
the  study  of  radiative  transfer  by  Ambarzumian  and 
Chandrasekhar,  [4].  More  generally,  it  utilizes  the 
"point— of— regeneration"  technique  of  the  type  used  by 
Bellman  and  Harris  in  the  study  of  branching  processes, 

[5l;  see  also  Harris,  [6]. 

Over  the  last  few  years,  frequently  in  collaboration 
with  Ueno  Ct] #  we  have  derived  a  large  number  of  functional 
equations  describing  a  variety  of  physical  processes,  and 
carried  out  some  large-scale  numerical  calculations 
([8],  [9]). 
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2 .  THE  MATRIX  RICCATI  EQUATION 

The  equation  we  wish  to  study  is 

(1)  R • (x)  =  B(x)  +  D(x)R(x)  +  R(x)D(x)  +  R(x)B(x)R(x) , 
R(0)  -  0, 

where  B,  D,  and  R  are  N  x  N  matrices  and  it  is 
assumed  that 

(2)  (a)  d1J(x)  >  0,  1  /  J, 

(b)  djj(x)  <  0, 

(c)  b^x)  2  °- 

Rather  than  tackle  this  equation  directly,  let  us  indicate 
its  physical  source,  and  then  show  how  the  simultaneous 
consideration  of  R(x)  and  two  related  functions  enables  us 
to  establish  existence  of  the  solution  of  (1)  for  all 
x  >  0  in  a  simple  and  painless  fashion. 

_ STEADY-STATE  NEUTRON  TRANSPORT  WITH  DISCRETE  ENERGY 

LEVELS 

Let  us  begin  by  describing  a  model  of  a  steady-atate 
transport  process  which  will  be  the  explicit  or  implicit 
source  of  many  of  the  analytical  ideas  we  shall  utilize 
in  what  follows. 

Consider  an  idealized  neutron-transport  process 
taking  place  in  a  one-dimensional,  homogeneous  isotropic 
rod  extending  along  an  axis  from  z  =  0  to  z  =  x.  We 


-3- 


suppose  initially  that  there  are  only  a  finite  number,  N, 
of  different  types  of  particles  moving  along  the  rod. 

These  possible  states  can  be  considered  to  be  energy 
levels,  labelled  i  =  1,  2,  •••,  N. 

It  is  assumed  that  when  a  particle  in  state  i 
traverses  a  segment  of  the  rod,  it  is  subject  to  inter¬ 
actions  with  the  substance  composing  the  rod.  These 
interactions  produce  two  possible  effects:  forward  or 
backward  scattering  into  any  of  the  N  possible  states, 
and  absorption.  However,  no  fission  occurs,  which  Is  to 
say,  there  Is  no  spontaneous  generation  of  new  particles. 

It  follows  that  the  total  number  of  particles  in  the 
process,  taking  account  of  those  absorbed  as  well  as  of 
those  scattered,  is  changed  only  by  addition  from  an 
external  source.  This  obvious  conservation  principle 
will  be  the  key  to  the  result  obtained  below  concerning 
the  existence  of  the  solutions  of  (2.1). 

In  this  paper,  we  shall  exclude  the  possibility  of 
collisions  or  Interactions  between  neutrons  themselves. 
This  will  permit  us  to  use  ordinary  differential  equations 
In  our  application  of  invariant  imbedding.  Subsequently, 
when  dealing  with  collisions,  we  shall  encounter  hyper¬ 
bolic  partial  differential  equations. 

Finally,  let  us  note  that  as  far  as  the  analysis  is 
concerned,  one-dimensional  transport  with  energy  levels 
is  equivalent  to  two-dimensional  transport  in  a  plane 
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parallel  slab  with  energy  and  angular  dependence.  We  are 
thus  treating  a  quite  general  transport  process  connected 
with  a  geometric  figure  such  as  sphere,  cylinder,  or 
plane-parallel  region. 


4_  l  _  .ANALYTIC  PRELIMINARIES 

Let  us  now  make  the  model  of  a  transport  process 
discussed  above  more  precise.  We  suppose  that  when  a 
particle  in  state  J  (J  a  1,2,**‘,N)  enters  the 
infinitesimal  segment  of  length  A  contained  in 
[x  +  A,x]  from  either  direction  (the  assumption  of 
isotropy),  the  following  events  take  place: 


(1)  (a)  The  expected  number  leaving  the  segment  in 

state  J,  moving  in  the  same  direction,  is 
1  +  +  0(A2). 

(b)  The  expected  number  leaving  the  segment  in 
state  1,  i  /  J,  moving  in  the  same  direction, 

p 

is  d^j(x)A  +  0(a).  (Forward  scattering.) 

(c)  The  expected  number  leaving  the  segment, 
moving  the  opposite  direction  in  state  i,  is 
blj(x)A  +  0(A2) .  (Back  scattering.) 

(d)  The  expected  number  absorbed  by  the  medium  is 
fjj(A)  +  0(A2). 


We  call  the  matrices  D(x)  =  (d^x)),  B(x) 
F(x)  -  (fii(x)6ij)»  the  forward  scattering,  back 


(t>ij(x)). 


scattering,  and  absorption  matrices,  respectively. 
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We  assume,  on  physical  grounds,  that 


(2)  (a)  d1J(x)  >0,  i  /  J, 

(b)  b1(J(x)  >  0, 

(c)  fu(x)  >  0, 


for  x  >  0.  The  basic  assumption  of  the  conservation  of 
matter  requires  that 


(3) 


djj(x) 


JlVX>  +  ifidU(x)  +  r33M  ■ 

1  Pi 

J  ■  1|2  j  * • * |Ni 


This  implies  that  <  0,  a  condition  required  to 

account  for  the  increments  to  other  states  and  for  the 
particles  absorbed. 

Introducing  the  matrix 


(4) 


\ 

1 

0 


°J 
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the  relations  of  (?)  can  be  written  in  the  simple  form 


(5)  M(B(x)  +  D(x)  +  F(x))  -  0,  x  >  0. 

This  is  the  fundamental  conservation  assumption 
which  will  yield  a  corresponding  conservation  relation 
for  the  matrix  functions  introduced  below. 


5.  REFLECTION,  TRANSMISSION  AND  LOSS  MATRICES 

Let  us  now  introduce  the  following  functions.  For 
i,J  =  1,2, ...,N,  let 

(l)  =  expected  flux  of  neutrons  in  state  i, 

reflected  from  a  rod  of  length  x, 
resulting  from  an  incident  flux  at 
x  of  unit  intensity  in  state  J; 

tij(x)  =  expected  flux  of  neutrons  in  state  i, 
transmitted  through  a  rod  of  length 
x  resulting  from  an  incident  flux  at 
x  of  unit  intensity  in  state  J; 

tlj(x)  =  expected  flux  of  neutrons  in  state  i, 
absorbed  within  a  rod  of  length  x, 
resulting  from  an  incident  flux  at 
x  of  unit  intensity  ir.  state  J. 

Schematically, 


(x) 


0 


Fig.  1 


When  we  say  a  rod  of  length  x,  we  mean  one  whose  ends 
are  respectively  at  the  fixed  position  0  and  the 


variable  position  x,  as  pictured  above.  As  a  conse¬ 
quence  of  the  assumption  made  above  concerning  no 
interaction  between  neutrons,  the  reflections,  trans¬ 
missions  and  absorptions  depend  linearly  upon  the  intensity 
of  the  incident  flux.  Hence  we  may  restrict  ourselves 
here  to  unit  incident  fluxes. 

Let  R(x)  =  (r1  j(x) ) ,  T(x)  =  (tjjCx)),  L(x)  »  t^U)) 

be  called  respectively  the  reflection,  transmission,  and 
absorption  matrices.  Using  invariant  imbedding  techniques, 
as  indicated  in  [l],  we  can  derive  differential  equations 
for  these  matrices.  For  the  sake  of  completeness,  let 
us  present  the  derivation  here. 

Consider  the  process  described  above  for  a  rod  of 
length  x  +  A  and  let  an  incident  flux  c  be  applied  at 
x  +  A.  Here  c  is  a  vector  flux  whose  J— th  component 
represents  the  intensity  of  the  incident  flux  in  state 
J.  By  virtue  of  the  assumptions  of  Sec.  4,  this  results 
in  a  flux  of  (I  +  Da)c  Incident  at  x,  a  flux  of  BcA 
reflected  from  x  +  A,  and  a  flux  of  FcA  absorbed,  all 

p 

to  terms  in  0(A  ).  The  flux  (I  +  Da)c  incident  at  x 
results  in  a  flux  R(x)(l  +  DA)c  reflected  at  x,  a 
flux  T(x)(I  +  Da)c  transmitted  through  the  rod,  and  a 
flux  L(x)(I  +  DA)c  absorbed  in  [x,0].  Here  B,  D, 
and  F  depend  upon  x,  as  indicated  above. 

The  flux  R(x)(l  +  DA)c  now  enters  the  segment 
[x,x  +  a]  and  results  in  further  interactions.  As  a 
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result  of  this,  we  have  the  additional  reflection 
(1  +  Da)R(x)(I  +  AD)c,  an  additional  absorption 
(af)R(x) ( I  +  DA)e  ■  AFR(x) c  +  0(a2),  and  a  quantity 
ABR(x)c  +  0(a2)  as  incident  flux  upon  [x,0].  This 
incident  flux  results  in  a  reflection  from  [x,0]  of 
R(x) ( ABR(x)c)  +  0(a2),  a  transmitted  flux  of 
T(x) { ABR(x)c)  +  0(a2),  and  a  loss  of  L(x) (ABR(x)c)  + 
0(A2).  The  flux  RBRA  +  0(A2)  through  [x,x  +  a] 
contributes  RBRA  +  0(a2)  to  the  total  reflected  flux. 

Schematically,  we  have 


( 1+DA ) R ( x ) ( I+DA ) c 


Fig.  2 


Adding  up  these  effects,  we  obtain  the  recurrence 
relations 

(2)  R(x+A)c  -  BcA  +  (I+DA)R(x)(I+Da)c 

+  R(x)BR(x)cA  +  0(A2), 

T(x+A)c  »  T(x)(l+DA)c  +  T(x)BR(x)cA  +  0(a2), 
L(x+A)c  -  PcA  +  L(x)(l  +  DA)c  +  FR(x)cA 
+  L(x)BR(x)cA  +  0(A2). 

Since  these  equations  hold  for  arbitrary  c,  we  can 
discard  c.  Expanding  the  left-hand  side,  and  passing  to 
the  limit  as  A  -*  0,  we  obtain  the  Riccati  differential 
equations 

(3)  R'(x)  -  B  +  DR(x)  +  R(x)D  +  R(x)BR(x), 

T'(x)  -  T(x)(D  +  BR(x) ) , 

L'(x)  -  L(x)(D  +  BR(x)  +  F( I  +  R(x) ) , 
with  the  physically  obvious  initial  conditions 

(4)  R(0)  -  0,  T(0)  -  I,  L(0)  -  0. 

Observe  that  of  the  three  functions,  it  is  only  the 
reflection  function  which  occurs  alone,  independently  of 
the  other  two.  The  remaining  two  functions  have  been 
deliberately  introduced  to  take  advantage  of  conservation 
properties . 


-10- 


6,  EXISTENCE  AND  UNIQUENESS  OF  SOLUTIONS 

The  conventional  existence  and  uniqueness  theory  of 
ordinary  differential  equations  establishes  the  existence 
and  uniqueness  of  a  solution  of  (50)  over  some  initial 
interval  [0,a].  Since  it  is  intuitively  clear  that  the 
reflection,  transmission,  and  loss  functions  must  exist 
for  all  x  >  0  (since  no  fission  is  allowed),  the 
question  arises  as  to  how  to  establish  this  analytically. 

For  the  case  of  constant  coefficients  (homogeneous  rod), 
we  can  reduce  the  equations  to  linear  equations  with 
constant  coefficients  and  use  the  explicit  solutions  to 
help  us.  For  inhomogeneous  equations,  this  approach  is 
more  difficult. 

To  obtain  an  equivalent  linear  equation,  let  us  first 
proceed  formally.  Consider  the  two  first  order  matrix 
equations 

(1)  X'  =  EX  +  FY, 

Y'  =  GX  +  HY, 

where  E,  F,  G,  and  H  can  be  dependent  on  x.  Consider  the 
matrix  Z  «  XY— ^ .  We  have 

Z1  =  (X1T1)1 

=  (EX  +  FY) Y— 1  -  X(Y-1(OX  +  HY)Y-1) 

=  EZ  +  F  -  ZGZ  -  ZH, 


(2) 
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an  equation  similar  to  that  satisfied  by  R(x).  The 
identification  is  complete  if  we  set 

(3)  E  -  D,  P  -  B,  H  -  -D,  G  *  -B, 
so  that  (l)  becomes 

(4)  X*  =  DX  +  BY, 
y»  =  — BX  -  DY. 

This  procedure  is  much  more  than  formal,  since  it 
turns  out  that  the  equations  of  (4)  are  the  transport 
equations  obtained  by  applying  the  usual  procedure  to 
the  study  of  the  fluxes  inside  the  rod. 

As  indicated  above,  it  is  not  a  trivial  matter  to 
study  (4)  when  B  and  D  are  variable  matrices. 

To  establish  nonlocal  existence  of  the  solutions  of 
(5*5) t  we  add  two  ingredients:  nonnegativity  of  the 
matrices  R(x),  T(x),  and  L(x),  and  the  conservation 
relation 


M(R(x)  +  T(x)  +  L(x) )  -  M. 

Both  of  these  conditions  are  intuitively  clear,  and,  as 
we  shall  see,  readily  established  rigorously.  Once  we 
have  done  this,  it  follows  that  R{x),  T(x),  and  L(x) 
are  uniformly  bounded  over  any  interval  of  existence. 

It  follows  that  the  solutions  can  be  continued  for  all 
x  >  0.  We  are  going  through  this  in  some  detail  since 


-12- 


the  same  line  of  reasoning  can  be  employed  for  many 
classes  of  functional  equations  arising  in  mathematical 
physics . 

7.  PROOF  OF  CONSERVATION  RELATION 

To  establish  the  conservation  relation  of  (6.1),  we 
consider  the  function 

(1)  Q(x)  -  M(R(x)  +  T(x)  +  L(x) ) 

and  differentiate  it  with  respect  to  x.  We  have 

(2)  Q'(x)  -  M(R'  +  T«  +  L») 

-  M(R  +  T  +  L)(TR  +  D)  +  M(DR  +  FR  +  B  +  F) 

-  Q(BR  -I-  D)  +  M(DR  +  PR  +  B  +  P), 

upon  using  the  equations  of  (5*3) • 

Considered  as  a  differential  equation  in  Q,  we 
observe  that  (2)  is  satisfied  by  Q(x)  ■  M,  since 

(3)  M(BR  +  D)  +  M(DR  +  PR  +  B  +  P) 

■  M{B  +  D  +  P)(R  +  I)  -  0, 

by  virtue  of  (4.5).  Since  Q(0)  =  M,  we  see  that 
Q(x)  «  M  within  the  interval  of  existence  of  R(x),  T(x), 
and  L(x).  This  argument  can  now  be  repeated  from 
interval  to  interval. 

It  is  remarkable  that  one  has  to  use  this  sophisti¬ 
cation  to  establish  a  relation  which  is  so  immediate 
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from  physical  considerations.  One  would  expect  in  place 
of  (2)  merely  the  relation  Q’  ■  0. 

8.  PROOP  OP  NONNEGATIVITY 

Local  existence  and  nonnegativity  of  solutions  can 
be  established  in  several  different  ways.  One  way  is  to 
convert  the  original  system  of  differential  equations 
into  a  set  of  integral  equations.  Let  us  begin  by 
writing  the  differential  equations  in  the  form 

(D  ^(e“DxR(x)e"Dx)  -  e_Dx[B  +  R(x)BR(x)  le""1*,  R(x0)  -  Rq 
^(T(x)e_Dx)  -  T(x)BR(x)e—Dx#  T(xq)  -  TQ, 

^(L(x)e_Dx)  -  [L(x)BR(x)  +  P  +  FR(x)]e_Dx,  L(xQ)  -  1^ 


Thus  an  appropriate  set  of  integral  equations  is 

D(x-x0)  D(x-xn) 


(2) 


R(x)  -  e 


R0e 


x  D(x— x^) 


D(x-x1) 


+  J  e  [B  +  R(x1)BR(x1) ]e  *  dXj 


*0 

-  01(R,T,L), 


D(x— xn)  „  x  D(x— x, ) 

T(x)  ■  TQe  0  +  J  T(x1)BR(x1)e  1  dXj 

*0 


02(R,T,L)f 
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L(x) 


D(x-x-) 

-  L0e 

x  D(x-x.  ) 

+  j  [L(x1)BR(x1)  +  F  +  PR(x1)]e  x  dx1 

x0 


™  ( Ri Tj  L)  • 


The  principal  result  we  wish  to  employ  to  establish 

Dx 

nonnegativity  is  that  d^j  >  0  implies  that  e  is  a 
nonnegative  matrix  for  x  >  Oj  see  [ 10] . 

Consider  the  space  S  of  triples  of  continuous 
matrix  functions  R(x),  T(x),  and  L(x)  defined  on 
xQ  <  x  <  xQ  +  a,  with  the  initial  values  R(xq)  -  RQ» 
T(xq)  -  Tq,  l(x0)  »  all  nonnegative  matrices, 

satisfying  the  constraints 


(3)  I  |R(x)  |  |  <  cx,  |  |T(x)  |  |  <  |  | L(x)  ||  <  Cj, 

where 


(4)  c1  >  Max  [||R0I|,  |  |  Tq  |  | , 

Consider  the  mapping  9  defined  on  S  by  means  of 
the  right-hand  sides  of  (2).  It  is  readily  seen  that  T 
is  a  contractive  mapping  of  S  into  itself,  provided 
that  a  is  sufficiently  small.  Thus,  by  virtue  of  the 
Cacciopoli  fixed— point  theorem,  0  has  a  unique  fixed 
point,  the  solution  of  (2). 

Alternatively,  we  can  construct  the  solutions  as 
the  limit  of  a  sequence  of  successive  approximations 
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given  by 

(5)  Rn+i  =  MVW'  n^°* 

Tn+1  =  02^Rn,Tn,Ln^ 

Ln+1  "  63^Rn,Tn,Ijn^ ' 

Applying  the  foregoing  result  with  Xq  »  0,  ci  >  N> 
where  N  is  the  dimension  of  the  system,  we  obtain  a 
solution  over  an  interval  0  <  x  <  a.  From  the  conserva¬ 
tion  relation  combined  with  the  nonnegativity  of  R(x), 
T(x),  L(x)  on  0  <  x  <  a,  it  follows  that  R(x),  T(x), 
L(x)  are  uniformly  bounded .  In  fact, 

(6)  I |R(x)| |,  I |T(x) | | ,  | | L(x) | |  <  N,  0  <  x  <  a. 

We  can  therefore  apply  the  result  with  xQ  »  a  and 
the  same  c^  as  before.  The  solution  can  thus  be 
continued  indefinitely. 

A  third  approach  starts  with  the  difference  equations 
obtained  from  (5.2)  by  neglecting  the  terms  which  are 

p 

0(A  ).  The  matrices  R(x),  T(x),  L(x)  are  defined  in 
this  way  for  x  -  0,  A,  2A,  ••*,  and  defined  by  means  of 
linear  interpolation  for  other  values  of  x.  Since 
I  +  DA  >  0  for  small  A,  we  see  that  the  matrices 
obtained  in  this  fashion  are  nonnegative  for  x  £  0. 

As  is  well  known,  these  functions  approach  the  solutions 
of  the  differential  equation  In  an  initial  interval 
0  £  x  £  b,  thus  once  again  establishing  nonnegativity. 
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9.  STATEMENT  OF  RESULT 

We  have  thus  established  the  following  result. 


Theorem.  If 


(1)  (a)  bjj(x)  >  °» 

(b)  djLJ(x)  >0,  1  /  J, 

(c)  fjj(x)  >  0» 


and 


M(B(x)  +  D(x)  +  P(x))  =  0 


for  x  >  0,  where  B  =  (b^fx)),  D  ■  (d^x)), 

P  =  ( f j j (x ) ^ »  and  M  "  then  the  equations 

(2)  R'(x)  -  B  +  DR(x)  +  R(x)D  +  R(x)BR(x),  R(0)  -  0, 

T'(x)  =  T(x)(D  +  BR ( x ) ) ,  T(0)  -  I, 

L»(x)  -  L(x)(D  +  BR ( x ) )  +  P(I  +  R(x) ) ,  L(0)  -  0, 


possess  a  unique  solution  for  x  >  0.  This  solution 
satisfies  the  conservation  relation 


(5)  M(R(x)  +  T(x)  +  L(x))  -  M, 
for  x  >  0. 

In  physical  terms,  this  means  that  the  reflection, 
transmission,  and  loss  matrices  are  defined  for  x  >  0, 
and  satisfy  the  equations  of  invariant  imbedding. 
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